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^ : ABSTRACT 

Results are presented of a harmonic analysis of the large scale cosmic-ray 
^ ■ anisotropy as observed by the Milagro observatory. We show a two-dimensional 

^ ■ display of the sidereal anisotropy projections in right ascension generated by 

^ ' the fitting of three harmonics to 18 separate declination bands. The Milagro 

observatory is a water Cherenkov detector located in the Jemez mountains near 

CN ■ Los Alamos, New Mexico. With a high duty cycle and large field-of-view, Milagro 
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is an excellent instrument for measuring this anisotropy with high sensitivity at 
TeV energies. The analysis is conducted using a seven year data sample consisting 
of more than 95 billion events, the largest such data set in existence. We observe 
an anisotropy with a magnitude around 0.1% for cosmic rays with a median 
energy of 6 TeV. The dominant feature is a deficit region of depth (2.49 ± 0.02 
stat. ± 0.09 sys.) xlO~'^ in the direction of the Galactic North Pole centered at 
189 degrees right ascension. We observe a steady increase in the magnitude of 
the signal over seven years. 

Subject headings: Milagro — Cosmic Rays — Anisotropy — Galactic Magnetic 
Fields — Heliosphere 



Introduction 



Observation of the sidereal large scale cosmic-ray (CR) anisotropy at energies of 1 - 100 
TeV is a useful tool in probing the magnetic field structure in our interstellar neighborhood 
as well as the distribution of sources. Cosmic-rays at these energies are almost entirely of 
Galactic origin and are expected to be nearly i sotropic due to interactions with the Galactic 
magnetic field (GMF) (IBerezinskii et al.lll990l ). The gyro radii of CRs at these energies in 
a GMF of about l/zG are from about 100 AU - 0.1 pc which is much smaller than the size 
of the Galaxy. This has the effect of trapping the CRs in the Galaxy for times on the order 
of 10^ years. Inhomogeneities in the GMF are interspersed randomly throughout the galaxy 
and in effect act as scattering centers for CRs. This randomizes their directions as they 
propagate leading to a high degree of isotropy on scales of a few hundred parsecs. 

Anisotropy can be induced through both large scale and local magnetic field configura- 
tions which cause deviations from the isotropic diffusion approximation. At lower energies 
of around a TeV, the heliosphere may be able to induce a C R excess in the direction of th e 
heliotail and also could modulate the overall CR anisotropy ( iNagashima et al.lll989l . Il998l ). 
At higher energies, the contributi on of discrete CR sources has been show n to be capable of 
creating a large scale anisotropy (jPtuskin et al.ll2006l : IStrong et al.l 120071 ). 



Diffusion of CRs out of the Galactic halo can also produce an anisotropy. Since the 
matter density is higher in the Galactic disk compared to that in the surrounding halo, the 
diffusion coefficient will generally be much higher in the halo. For this reason, CRs pro- 
duced in the Galactic disk will tend to diffuse out into the halo creating an anisotropy in 
the direction perpendicular to the disk. Predictions of the anisotropy have been made using 
values of the diffusion coefficient inferred using a given propagation model and observational 
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data. This predicted anisotropy can take on values of 10~^ to 10~^ depending on the prop- 
agation model used (IBerezinskii et al.lllQQOl ). Given this correlation between the anisotropy 
and diffusion coefficient, knowledge of the large scale CR anisotropy can be used to constrain 
diffusion models. 



In addition to the above effects, Compton and Getting introduced a theory (ICompton fc Getting 



19351 ) of CR anisotropy which predicts a dipole effect due to the motion of an observer with 
respect to an isotropic CR plasma rest frame. The anisotropy arising from the Earth's mo- 
tion around the sun is calculated to be on the order of 10"'' and wou ld appear in universa l 
time. This effect has bee n observed by numerous experiments (e.g. see lAglietta et al.l (119961 ) . 
Amenomori et al.l (120041 )). There is also a possible sidereal time anisotropy coming from the 
motion of the Solar System through the Galaxy. This effect is more difficult to predict given 
that the isotropic CR rest frame is not known. Using the assumption that the CR rest frame 
does not co-rotate with the Galax y leads to a predicted an isotropy on the order of ~ 0.1%. 
This effect has not been observed (lAmenomori et al.ll2006l ). 



There have been numerous observations of the large-scale sidereal anisotropy in the range 
of energies 10^^ — 10^^ eV. In this paper we use large-scale anisotropy to refer to features 
of the anisotropy in the sky with an extent greater than ~ 40° in right ascension. Most 
of these observations are examined by fitting harmonics to the distribution of cosmic-ray 
events in the right ascension direction. The main features of the anisotropy, the amplitude 
and phase, are fairly constant over the energ y range mentioned ab ove (3 — 10 x 10~^ for 
the amplitude, and — 4 hr for the phase, see iGuillian et al.l (120071 ) for a compilation). At 
TeV energies, the focus of this paper, the anisotropy is known to have a value on the order 
of 10~^ with a deficit region, sometimes called the "loss-cone", around 200° right ascension. 



and an excess, or "tail-in" region, around 75° dNagashima et al.lll989l . ^9981 ). In this paper 



we present the results of a harmonic analysis of the large scale cosmic-ray anisotropy in the 
northern hemisphere as observed by the Milagro experiment. 



2. The Milagro Observatory 



The Milagro observatory (lAtkins et al.l |2004| ) is a water Cherenkov detector which is 
used to monitor extensive air showers produced by TeV gamma-rays and hadrons hitting the 
Earth's atmosphere. Milagro is located in New Mexico at a latitude of 35.88°, a longitude 
of 106.68°, and an altitude of 2630 m above sea level, possessing a large field of view of ~2sr 
and a high duty factor of > 90%. The detector is composed of an 80m x 60m x 8m pond filled 
with ~ 23 million liters of purified water and protected by a light-tight cover. The central 
pond is instrumented with two layers of PMTs: a top layer with 450 PMTs under 1.4m of 
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water which detects Cherenkov hght from air shower electrons, electrons Compton scattered 
by gamma-rays, and gamma-rays that have converted to electron-positron pairs in the water; 
and a bottom layer with 273 PMTs 6m under the surface used for gamma-hadron separation 
(not used in this analysis). The direction of an air shower is reconstructed using the relative 
timing of the PMTs hit in the top layer of the pond with an angular resolution of < 1°. 

This pond is surrounded by a 200m x 200m array of 175 "outrigger" tanks. Each "out- 
rigger" is a cylindrical, polyethylene tank with a diameter of 2.4 m and a height of 1 m. 
The outrigger tanks contain ~ 4000 liters of water and are instrumented with a single down- 
ward facing PMT located at the top of the tank. The inside of each tank is also lined with 
Tyvek (a white, reflective material) to increase the light collection capability of the PMT. 
The outrigger array, completed in 2003, is used to improve the angular resolution and in the 
estimation of primary particle energy. 

The data used in this analysis were collected by Milagro from July 2000 through July 
2007. The hardware trigger was essentially a voter coincidence, requiring a multiplicity of 
about 75 PMTs. During this time the trigger multiplicity was tuned to maintain an average 
trigger rate of ~1700 events per second (to within 10 %), the majority of which are due 
to hadronic showers. After event reconstruction, accepted events are required to have used 
at least 50 PMTs in the angular fit {Nfu) and have a zenith arrival angle of < 50°. The 
Nfit > 50 requirement removes effects due to fluctuations in the trigger threshold as events 
passing this cut generally have a higher PMT multiplicity than the hardware trigger requires. 
Although this makes the data more stable, it also reduces the effective trigger rate to ~ 500 
Hz. After these cuts the data set consists of 9.59 x 10^° cosmic-ray events with a median 
energy of 6 TeV as determined from the measured primary cosmic-ray spectra and simulated 
detector efficiency. The median energy varies slowly with zenith angle; it is 4 TeV for zenith 
angles from 5° to 10°, and 7 TeV for zenith angles from 35° to 40°. 

Detector response to primaries is determined with a Monte Carlo simulation. This 
Monte Carlo has two main components: the simulation of an air shower in the atmosphere 
and the simulation of the showe r through the M ilagro detector. These simulations are 



carried out using the CORSIKA Jfleck et al.lll998h and GEANT4 jAgostinelh et al.ll2003h 



packages respectively. Cosmic-ray nuclei ranging in mass from the dominant protons to 
those as hea vy as iron are sirn ulated according to their spectra as measured by the ATIC 



experiment (jPanov et al.ll2006l ). Because of the large shower to shower variations, Milagro 
cannot measure the energy of individual primaries generating these air showers without 
careful data selection criteria. Measurement of energy correlated parameters does allow 
some determination of primary energy spectra without severely restricting the data. The 
energy analysis of the large-scale anisotropy will be the subject of a future publication. 
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3. Data Analysis 

3.1. Forward-Backward Asymmetry and Harmonic Method 

We have designed an analysis method to search for fractional deviations from isotropy 
down to levels of about 10~^, below the level of predictions of the Compton-Getting effect 
and the magnitude of previous observations. Such sensitivity cannot be achieved by looking 
at raw event rates as a function of direction, because the air shower development essentially 
makes the atmosphere part of the detector, and weather effects lead to typical overall trig- 
ger rate variations of ±15%. Thus, counting raw event rates would give random apparent 
anisotropics completely overshadowing a true signal. Instead, we employ the technique of 
forward-backward asymmetry (FB) using the number of events {Np and Nb) collected in 
some small time interval in two "telescopes" , i.e. pixels of small and equal solid angle, at the 
same forward and backward angle as shown in Figure [H 



FB = {Nf - Nb)/{Nf + Nb) (1) 

The expression of FB is manifestly independent of overall detector rate, as can be seen 
from the invariance of FB under the substitution Np, Nb cNp, cNb, where c is a constant. 

The analysis method utilizes the rotation of the Earth to search for a coherent modula- 
tion of FB during a 24 hour day. The FB modulation does not directly give the anisotropy 
of the sky, but is a function of it. As a region of the sky with an excess cosmic-ray flux 
relative to the average (a positive anisotropy region) is swept into the forward "telescope" 
FB becomes more positive; when the same excess is swept into the backward "telescope", 
FB becomes negative. One can think of FB as a coarse "derivative" of the actual anisotropy 
of the sky. Thus the FB modulation is a tool to obtain the quantity of interest, the fractional 
anisotropy (i.e. deviation from uniformity) of the sky, as detailed below. 



The FB method is commonly used in particle physics (e.g. see iBand et al.l (Il989l )) and 
is also closely related to the Ea st- West technique , classically used in other anisotropy mea- 
surements (for an example see Nagashima et al. ( 19891 )). in which the ratio examined for 



modulation is EW = {Ne — Nw)/Ntot, formed with the rate integrated over the whole East- 
ern and Western sky. Again this is a rate independent quantity which is used to remove 
random fluctuations in overall detector rate, a common requirement in all experiments using 
the atmosphere as part of their detector. The EW method can be thought of as a single 
integrated measurement, compared to our multiple localized and independent measurements 
of FB which are used to improve the statistical power of this analysis. 



Random daily weather-induced or instrumental variations of the anisotropy itself (and 
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therefore of FB) are averaged out by summing (i.e. averaging) many full days, each sweeping 
a full circle in the sky. The random anisotropies decrease as N^ays while a coherent signal 
is not reduced. Finally, because the method uses the time modulation of FB rather than 
FB itself, it is impervious to the inherent geometrical anisotropy present in the detector 
acceptance seen in the azimuthal distribution of Figure [2l which was observed to be stable 
during the lifetime of the experiment. We note that only instabilities of this geometrical 
anisotropy on the scale of a single day would affect the measurement, while a slow evolution 
would not. 

Because this method measures the modulation in the direction of the Earth's rotation, it 
yields no information about the modulation in the declination (dec.) direction. The results 
will therefore be for the projection of the anisotropy in the right ascension (r.a.) direction 
rather than the full 2-D anisotropy of the sky. Such a projection can be created for any visible 
dec. band. Since each dec. band sweeps around a full circle of the sky independently and 
contains statistically independent data, we choose to do separate analyses for each 5° dec. 
band of our data, considering each as a separate observation. Any theoretical description 
of the true 2-D anisotropy can be confronted with, or constrained by, our data (given in 
Section 4, Table 1) by projecting the anisotropy along r.a. in our dec. bands. In the rest of 
this paper we will always use the word anisotropy to mean this projected anisotropy. 

To do the harmonic fit we make the assumption that the large scale anisotropy in any 
given dec. band can be modeled by a Fourier series and that it is a small modulation of 
a nearly isotropic signal. Three harmonics (the fundamental and the next two) have been 
found to be sufficient for this method (see section 4.1). This allows us to see large scale 
effects having a width in r.a. of greater than ~ 40°. In this harmonic model the normalized 
rate in a given direction of the sky, in celestial equatorial coordinates with declination = 5 
and right ascension = 9 [oi the equivalent coordinate in the right ascension direction when 
using a different coordinate system) is: 



where for a given 5 the average in the denominator of the left side is over 6', and A^iO) is 
the fractional anisotropy in the right ascension direction. This model of the sky is complete 
once the six Fourier coefficients (three amplitudes and three phases) on the right side of the 
equation are known. The data analysis thus consists of determining these Fourier coefficients. 

It is this harmonic model that allows us to connect the quantity of interest, namely 
the fractional anisotropy in the sky As{9), and our chosen tool, the measurement of the 




n=l 
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modulation of the forward-backward asymmetry. This connection is detailed exphcitly below. 

The quantity Ng^^sl^) is the number of cosmic-ray events collected during a particular 
time interval (characterized by the angle 6q), in an angular bin at a given declination (6) 
and local hour angle (^) characterizing the symmetric forward (+0 ^ind backward (— 
inclination of the "telescope" of Fig. [H In this notation the forward-backward asymmetry 
of Eq. 1 becomes: 



(3) 



The bin counts are computed for 1/ 2-hour intervals in sid ereal time (ST), universal 
time (UT), and anti-sidereal time (AST) (IFarley fc Storeyill954j ). These 1/2-hour intervals 
are parameterized by an angle Oq which specifies the relative advance of the local meridian 
through the sky for the three different time scales: 



= 3.75° + 7.5° X 1ST 



00 = 3.75° + 7.5° X lUT 
Oo = 3.75° + 7.5° X I AST 



1ST, JUT, and lAST are integers, from zero to 47, denoting half hour intervals of sidereal 
time, UT, and anti-sidereal time (defined by flipping the sign of the transformation from 
universal time to sidereal time). In the above equations, the constants convert an integration 
time interval (1/2 hour) into degrees (7.5°) and a starting angle at the center of that interval. 

In practice the cosmic-ray events are recorded in histograms with 5° x 5° bins according 
to their arrival direction from —10° to 80° in declination and —50° to +50° in hour angle. 
The events are collected over a 30 "minute" period giving us 48 half hour histograms per day, 
where "minute" is defined in the three following time frames: sidereal (366.25 days/year), 
universal (365.25 days/year) and anti-sidereal (364.25 days/year). The coordinate system 
in the sidereal time frame is "sky-fixed", in the universal time frame it is "sun-fixed", and 
in the anti-sidereal frame it is "non-physical" . Histograms for a given half hour period are 
accumulated over any number of days to build an average set which corresponds to a chosen 
time period. A representative example is shown in Figure [3l These 48 half hour histograms 
for the time frame, period, and dec. band of interest are then analyzed using the method of 
forward-backward asymmetry described above combined with a harmonic fit. 

Recognizing that 6 = 6q ± ^, substituting Eq. [2] for the harmonic sky model into Eq. [3l 
and then utilizing 7 <C 1 and the appropriate trigonometric identities, we get: 
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3 



FBs{9o, ~ X] -7n,<5 sin«) sm(n(6'o - 0n,5)) 



(4) 



n=l 



Equation (jl]) is the prediction of the harmonic model of the sky, in terms of its Fourier 
coefficients, for the forward-backward asymmetry, to be fit with our experimental data. 

The next step is preparing the experimental data for the fit. For a fixed "telescope" 
angle ^ and at a given 6 we calculate FB and its statistical error by using the relevant pair 
of histogram bins in each of our 48 saved histograms of hour angle vs. dec. (see Fig. [3] for 
an example). Since we are only interested in the modulation of FB with (i-e. with time), 
the average over all Oq is subtracted from the calculated FB. This step separates FB due 
to the inherent geometric anisotropy of the detector from the time modulated FB due to 
the sky anisotropy. These values are assembled into a 48 point, one-dimensional histogram 
of FB and its error vs. for this particular choice of ^ and 6. Equation H] could be fit to 
this histogram to obtain the Fourier coefficients. But, this procedure can be done for several 
choices of C, that we call ^j, each is sampling the same sky, yielding a new FB vs. 6q graph 
to fit. By choosing a set of (where C,i ranges from 2.5° to a dec. dependent maximum of up 
to 47.5° in 5° steps), centered on each of the 5° x 5° pixels (see Fig. [T]&[3]) and calculating the 
FB (Eq. 3) with the contents of the corresponding pixels, we obtain up to 480 (48 values of 
9o times up to 10 values of ^j) statistically independent measurements per dec. band. Every 
available pixel (histogram bin) pair in the 48 hour angle vs. dec. histograms (Fig. [3]) is used 
once and only once. These one-dimensional histograms are assembled into a two-dimensional 
histogram of FB vs. and ^j. This histogram contains all our information about FB and 
its statistical error. The fit of Eq. IHto this 2-D histogram (a simultaneous fit over and ^i) 
gives the experimental values and errors of the Fourier coefficients for the harmonic model 
of the anisotropy in a dec. band (Eq. [2]). An example of the 2-D histogram that is being 
fit, and the result of the fit in 6q for a single slice in are shown in Fig. IH^a) and (b). The 
coherent signal in FB is clearly seen in both. Figure ID^b) also shows the resulting measured 
anisotropy reconstructed using the six fitted Fourier coefficients in this example. 

The determination of the statistical errors on the Fourier coefficients is straight-forward. 
The statistical errors on the number of events {Np and Nb) are simply their square roots 
and these errors are propagated to the error on FB: 



to a superb approximation. Since each FB bin in the 2-D histogram is statistically 
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independent of the others, the fit is that of a 2-D experimental histogram with independent 
errors to the parameters of the theoretical function of Eq. (jlj). This fit (using MINUIT 



( lJamesll2000l )) propagates the errors on the experimental data into errors on the parameters 



of the theory, namely the six Fourier coefficients of the anisotropy. 

It is noted that the simultaneous fit of the different C,i bands is in effect an averaging 
over which means that we are also averaging over any difference of energy in the data due 
to the dependence of atmospheric depth on ^. This averaging is well justified after the fact 
by an examination of the ^ dependence discussed in Section 14.31 

Finally, we emphasize that the role of the forward-backward asymmetry is over once 
the fit is done. The coefficients obtained are for the anisotropy of the sky, defined as the 
fractional difference of the rate at a given r.a. from its average over all r.a., for each dec. 
band. In Section UTTl we tabulate the fit results for these 18 dec. bands, as well as graphically 
display their behavior. 



3.2. Monte Carlo Tests of the Analysis Method 

Tests of the analysis method were conducted using a Monte Carlo simulation that takes 
as input a 2-D map containing any desired anisotropy in universal time (UT) or sidereal 
time (ST) and then outputs events which can then be analyzed using the method discussed 
in the previous section. 

We describe one important case in detail: a fixed UT signal will average to zero in 
ST or AST, but a seasonally modulated day-night anisotropy in UT will generate a false 
"sideband" signal in ST and AST. We use a time varying input in UT, where the magnitude 
of the UT signal varies sinusoidally between zero and its maximum with a period of one 
year, simulating an extreme seasonal variation. Figure [5] shows the results of an all-dec-band 
analysis in three time frames for this test. The ST and AST frames, which have no input 
signal in the MC, show a clear signal induced by the time varying UT input. In this example 
the induced signal is at about half the amplitude of the UT time variation. The phases of 
the signals in ST and AST are different but their magnitudes are the same. Using this result 
we can estimate the systematic errors in the sidereal signal by examining the anti-sidereal 
signal. Since there are no physical processes which occur in the anti-sidereal time frame the 
signal should be zero when data sets of an integral number of years are used as seen from 
the Monte Carlo checks above. Therefore a signal present in the non-physical AST frame is 
interpreted as being induced by variations in time of the UT signal and such a signal will 
also appear in the ST frame. 
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Two other sets of MC tests were successfully passed by the FB method. The first set 
tests the stability of the analysis method by using various simple signals (e.g. isotropic sky or 
the result of the analysis of real data) as input in ST with no input in UT or anti-sidereal time. 
Analysis of the output of these tests is consistent with the input within statistical errors. 
The UT and anti-sidereal time (AST) maps in these tests are consistent with isotropic as 
long as an integral number of years is used. This is to be expected since a static point in 
ST moves across the UT frame and returns to its original position after one year. In AST 
this point will move across two times in one year. Therefore when the declination band is 
normalized a sidereal signal will average to zero in both UT and AST. This cancellation is 
of course also true for the sidereal sky with a signal fixed in UT (e.g. day, night effects). 

The second set of tests shows that no false harmonic features are induced from isolated 
structures (e.g. a squar e well deficit cente r ed around ^ 185° or th e known Galactic equator 
excess of gamma-rays (lAtkins et al.ll2005l : lAbdo et al.l 120071 . l2008l )). These tests show that 
no extra signals are induced except in the cases of a sharply discontinuous input and a time 
varying UT input. The discontinuous input induces features due to the use of only three 
harmonics as is expected from Fourier theory. 



4. Results 

4.1. Multi-dec-band Results and the Sidereal Sky Map 

Table I is a summary of the harmonic fit parameters for the sidereal sky fractional 
anisotropy obtained in each of 18 individual dec bands. The for the fit and the sample 
size in each dec. band is also listed in this table. Three harmonics were chosen as they give 
a better x^/ndf than one or two harmonics whereas there was little improvement in using 
four. Fig. E] displays the anisotropy profiles in RA for eight adjacent dec. bands in the table. 
As seen from the definition of the anisotropy in Eq. [21 the profiles show the deviation from 
their average over all r.a., so that the area below and above the reference level of zero is the 
same: the existence of a deficit in some region implies an excess elsewhere. 

Figure [7] assembles all 18 anisotropy profiles in their respective dec. positions into a 
2-dimensional color display of the local anisotropy, with red representing an excess and blue 
representing a deficit. Care must be taken in interpreting Fig. [71 As explained before, our 
measurement is in the r.a. direction only, so that this picture does not purport to be the 
complete anisotropy of the celestial sphere; nor can the complete anisotropy be inferred from 
our results. Fig. [Tjis a two-dimensional display of one- dimensional information, the fractional 
anisotropy in r.a. What can be fairly compared across the dec bands are the shape and 
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strength of the fractional anisotropy in the r.a. direction; there is no information on the r.a. 
averaged cosmic ray rate difference from one dec. band to another, i.e. the anisotropy in the 
dec. direction. 

The data and their analysis are completely independent for each dec. band, but a striking 
commonality of the anisotropy behavior is seen in both Figs. [6] and [3 Adjacent dec. bands 
have very similar profiles, in shape and in strength. The conclusion is that there is indeed a 
coherent anisotropy signal over a large swath of the sky. The dominant feature is a prominent 
valley or deficit region extending from 150° to 225° in r.a., and clearly visible in all dec. bands 
between —10° and 45°. The decrease of the depth of this valley, towards large values of dec, 
seen in the color picture of Fig. [TJ is explained at least in part by the fact that the circle 
in the celestial sphere generated by the earth's rotation gets smaller and smaller, shrinking 
to a single point as dec. approaches 90°. The r.a. position of the valley minimum appears 
dec. independent in Figs. |6] and [71 Figure [8] plots the amplitude and phase of the lowest 
harmonic versus dec, confirming both these trends; the amplitude is seen to decrease as dec. 
gets larger exhibiting the expected cos(dec) dependence, while the phase is quite stable for 
dec. below 60°, where the amplitude is large enough for the phase to be well defined. 



4.2. Single-dec-band Results for ST, UT, and AST Time Frames 

As a complement to the multi-dec-band analysis of section 4.1 we do a single-band 
analysis by projecting the hour angle vs. dec. histograms onto the horizontal axis of Fig. [3l 
thus integrating all dec. bands, then performing the same anisotropy analysis on this single 
band, to obtain the overall r.a. behavior. Examining the resulting single band anisotropy, 
with only 6 fourier coefficients, gives a convenient and economical way to compare and 
contrast the behavior in all three coordinate systems corresponding to sidereal, universal 
and anti-sidereal time. Table 2 gives the fitted harmonic fit coefficients for ST, UT, and 
AST, and Fig. [9] displays the resulting profiles for all three coordinate systems. 

The sidereal-frame profile of Fig. M clearly reproduces the dominant feature seen in the 
multiband analysis, the valley centered at 189° in r.a. The single band valley depth (SBVD) 
is an alternate measure of the strength of the sidereal anisotropy and is found to be SBVD 
= (2.49 ± 0.02 stat. ± 0.09 sys.) xlO~'^, a 22.6a signal. The method of systematic error 
estimation will be discussed in Section 14.31 

The Compton Getting effect due to the earth's motion around the sun gives a predicted 
anisotropy in universal-time frame which is a dipole (the lowest harmonic), plotted in Fig. [9] 
together with our experimental result. Because the predicted signal amplitude is a factor of 
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five smaller than the observed sidereal time signal, and is correspondingly more susceptible to 
systematic distortions, this is a stringent test of the FB method. The experimental anisotropy 
profile is indeed dominated by the lowest harmonic as expected, has the predicted amplitude 
and reasonably close, though not in agreement (within errors) with the predicted phase. A 
detailed analysis of the UT measurement, including systematic errors, is deferred to a future 
publication. We consider our observation of the Earth-motion Compton-Getting anisotropy 
signal in UT a powerful check of the analysis method, and thus a verification of the reliability 
of the result for the sidereal sky. 

Since the anti-sidereal reference frame of Fig. [9] is non-physical, we expect no real signal 
there, while systematic effect distortions will show up in this frame also. Indeed this signal, 
plotted in Fig. [HI is seen to be a factor of 3 smaller than the UT signal and thus a factor of 
15 smaller than the ST signal. The observed AST signal will be used in the estimation of 
systematic errors in the next section. 

4.3. Stability of Data and Systematic Errors 

As a verification of the threshold independence and overall detector rate independence 
of the method and results, we raise the PMT multiplicity (the number of PMTs hit during 
an event in the air shower layer of the detector) requirement in the trigger from its hardware 
value of around 75 in several steps up to a value of about 400. Figure [10] shows the num- 
ber of events and the fractional anisotropy valley depth SBVD as a function of the PMT 
multiplicity. The valley depth is seen to be stable over a factor of five in trigger threshold 
corresponding to a factor of ~ 25 in trigger rate. This is a direct confirmation that the 
result is indeed impervious to the detector rate and threshold, as expected from the rate 
cancellation of the FB method, explained in Section 13.11 

In section 3.1 it is mentioned that the analysis method averages over different values of 
1^, from 2.5° to 47.5°. Figure [TT] compares the single-band profiles when the full data set is 
split in two and analyzed separately, where one set contains data with ^ ranging from 0° to 
25° and the other 25° to 50°. The good agreement of these two profiles demonstrates that 
the procedure of simultaneously fitting over all ^ is well justified. 

To check for the possibility of a dependence on seasons, we break the data into three 
seasonal sets, defined as: Apr. -Jul., Jul. -Nov., and Nov.-Apr. corresponding to average lo- 
cal weather periods of: warm with low precipitation, high precipitation, and ice and snow 
respectively. No significant changes are seen in Fig. [11] which shows the three single-band 
profiles superimposed. 
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To test whether the observed ST signal is a sidereal phenomena and not due to universal 
time effects, such as weather, we break the data into two month periods. On these short 
time scales there will not be a cancellation between the ST, UT and AST time frames since 
a fixed position in ST corresponds to a nearly constant position in UT. Fig. [T2] shows the 
position of the observed minima in UT for each two month period over seven years. The 
solid line shows the expected position of a fixed sidereal minimum located at 189° r.a. If 
the minimum observed at 189° r.a. is a dominant feature of the sky, then the minimum 
in UT should follow the sawtooth pattern of the solid lines in the figure. If, on the other 
hand, there were dominant features in the UT sky, the positions of the minima should not 
correlate at all with this sawtooth pattern. The observed correlation between the observed 
and predicted positions of the minima in Figure [12] is a strong indication that the ST signal 
is the dominant feature in UT on this time scale. The same conclusion is reached when the 
corresponding analysis is carried out in the AST time frame where the sawtooth pattern has 
double the frequency. 

To estimate systematic errors in ST or UT, for data sets of a full year (or several full 
years) of data, we use the anisotropy profiles in AST. The phase in the nonphysical anti- 
sidereal frame sweeps continuously between 0° and 360° in ST or UT in the course of a 
year, so that any coherent signal is lost because its phase is in effect randomized by this 
sweep. Since the physical signals have been washed out, what is left are any systematic 
distortions. We illustrate the logic and procedure for the case of the AST profile seen in 
Fig. M for the single band analysis. While the deviation of this profile from zero gives the 
possible magnitudes of the distortion of the valley depth, the phase relative to the signal 
in ST is unknown. If a deficit region of distortion in the AST profile were to coincide with 
the valley it would increase the valley depth SBVD, an excess would decrease it, and a zero 
crossing would leave it unchanged. We account for all these possibilities by projecting the 
AST curve of Fig. M onto the vertical axis, to obtain the histogram shown in Fig. [T3] of all 
the possible distortions of SBVD that can result from this profile. We then take the RMS 
deviation of this histogram, aRMS = 0.09 x 10^^ as the systematic error on the valley depth 
SBVD in the ST analysis. 



4.4. Time Evolution of the Sidereal Anisotropy 

To examine any time evolution of the anisotropy the data were split into seven year-long 
sets from July 2000 to July 2007. Figure [H] shows the three amplitudes and three phases 
obtained by fitting data from each yearly set using the single band analysis (over all dec). 
There is a stability in phase of all three harmonics over these seven years and in amplitude 
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of the two higher harmonics. The fundamental harmonic, however, shows a clear increase in 
amplitude with time. 

In order to quantify this increase in a precise manner, we look at the time evolution 
of the single band valley depth (SBVD) in the all dec. analysis, defined in Sec. 14.21 as a 
single-parameter measure of the strength of the anisotropy. Using this parameter allows us 
to include all three harmonics as well as utilize the systematic error estimation procedure 
outlined in the previous section. First we confirm that the valley is stable in position over 
time. The position in r.a. of the minima over the course of seven years shows very little 
variation as is expected from the stability of the harmonic phases seen in Fig. [TH Fitting 
the position of the seven yearly minima to a constant yields a value of 189° ± 1° r.a. with a 
'X^ /ndf = 4.5/6. This lack of change in position over time is what one would expect from an 
actual sidereal signal. 

Figure [151 SBVD vs. year, shows that there is strong evidence of a strengthening of the 
valley depth over the seven year span of this data set. 

To test the robustness of this time dependence a number of checks were done. As a test 
that is completely different from the insensitivity of anisotropy strength to trigger thresholds 
(described in Section H73] and Fig. [10]), we have done a direct check of whether the time- 
dependence of SBVD itself is threshold dependent. For several raised multiplicity thresholds 
between 90 and 280 PMTs hit in the top layer of the pond, the same time dependence is 
seen; the yearly trend does not disappear. 

To see that this is a sidereal effect and not a detector effect we look at the yearly 
time evolution for the universal time and anti-sidereal time signals. Figure [16] shows the 
amplitudes of the three fit parameters for the single band analysis (all dec.) in both UT and 
AST; the earth motion Compton- Getting effect in UT should have no time dependence of the 
amplitude. The amplitudes of the harmonics in UT are constant over this seven year data 
set, within the errors, as well as their phases (not shown in the figure). With respect to the 
amplitudes of the harmonics in AST, these appear to be significantly larger in some years, 
but even the largest amplitudes are 5 to 10 times less than those seen in ST. From these 
tests it thus appears that time dependent detector effects cannot account for the observed 
strong time dependence of the sidereal anisotropy. 

5. Conclusions 

Previous experiments such as the Tibet Air Shower Array, with a modal energy of 3 TeV, 
and Super-Kamiokande-I, with a median energy of 10 TeV, have identified two coincident 



- 15 - 



regions of interest in their sidereal observations: an excess located at ~ 75° r.a. or "tail- 
in" a nisotropy, and a de ficit at ~ 200° r.a. or "loss-cone" anisotropy (lAmenomori et al. 



20061 : iGuillian et al.l 120071 ). Both of these regions are consistent with Milagro observations. 
The "loss-cone" coincides with the deep central-deficit region seen in this analysis while the 
narrow "tail-in" excess is clearly observed in anoth er Milagro analysis which is sensitive to 



features with extent smaller than ~ 30° in r.a. (see lAbdo et al.l (120081 ) region A). 



The strengthening of the signal in the central-deficit region over time is a result unique 
to this analysis. Tibet found no evidence of time variation comparing data split into two 
five-year periods, 1997-2001 and 2001-2005. However, only the second of these time periods 
overlaps with our data set for which the average value we observe in the deficit region agrees 
well with their measured deficit. 

The anisotropy observed in the galactic cosmic rays could arise from a number of possible 
effects. The Compton-Getting effect (CG) predicts that due to the motion of the solar system 
around the galactic center through the rest frame of the cosmic-ray plasma an anisotropy 
is induced with the form of a dipole with a maximum in the direction of motion. For no 
co-rotation of the cosmic ray plasma with the Galaxy, the magnitude of the anisotropy is 
calculated to be 0.35% given our speed of ~ 220kms~^ , while at the other extreme of full 
co-rotat ion, the anisotropy wou ld be zero. No evidence of a Galactic CG anisotropy was 
seen in lAmenomori et al.l (120061 ). For no co-rotation, the dipole should have a maximum 
at r.a. = 315° and dec. = 48°, and a minimum at r.a. = 135° and dec. = —48°. Since our 
analysis method yields only a projection of the anisotropy the observed CG effect will be 
slightly different from the true effect. The CG effect we expect to see in this analysis was 
determined from Monte Carlo simulation and found to be a dipole with a maximum of 0.14% 
at 315° r.a. for the declination range between 50° to 60°. This range was considered to try to 
reduce effects from the central-deficit region. For the actual data, fitting a single harmonic to 
the projection in r.a. corresponding to declinations 50° to 60° gives a /d.o.f. = 11505/998 
which is clearly a poor fit. Although this suggests that the observed anisotropy is not 
dominated by the galactic Compton-Getting effect, its contribution to the anisotropy cannot 
be ruled out. 

In addition to the Compton-Getting effect there is expected to be an anisotropy stem- 
ming from the diffusion of cosmic-rays in the interstellar medium. At high energies the 
main effects are expected to be mainly due to the distribution of discrete CR sources and 
the structure of the galactic magnetic field in the neighborhood of the solar system. One 
study conducted consists of a simple diffusion mo del assuming increased production in th e 
galactic disk due to supernova remnants (SNR) (jPtuskin et al.l l2006l : IStrong et al.l 120071 ). 
Also considered was the diffusion of CR out of the galactic halo. This is attractive since 
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we have a deficit region in the general direction of the North Galactic Pole which could be 
due in part to this diffusion. The main contribution of CR from SNR was considered for 
sources with distances from Earth of < 1 kpc and ages < 0.05 Myr. Calculations performed 
(jPtuskin et al.ll2006l : IStrong et al.ll2007l ) using these sources and taking into account CR re- 
acceleration as well as diffusion out of the galaxy gives an anisotropy about 3 times greater 
than observed, with the main source of the anisotropy due to the Vela SNR located at 128° 
r.a. and —45.75° dec. This model only predicts the magnitude of the expected anisotropy, 
not its exact phase. It also is remarked that this is a simplified model which assumes an 
isotropic diffusion tensor which is not explicitly known to be true at these energies. 

At energi es of ~ 1 TeV the heliospher e is believed to have an influence on the distri- 
bution of CR jNagashima et al.lll989l . llQQsh . One possible reason for the modulation of the 
anisotropy on the observed time scale could be due to variations in the heliosphere since 
we know that it changes in relation to solar output. It is noted that our data begins at the 
solar maximum and ends near the solar minimum. A recent derivation of the diffusion tensor 
contains a new component due to perpendicular spatial diffusion which is expected to be 
an important factor in understanding the an isotropy due to the Gala ctic disk as well as the 
modulation of CR in the outer heliosphere (jSchlickeiser et al.l 120071 ). Finding a consistent 
explanation of the observed anisotropy and especially its time evolution will be a challenge. 
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Hour Angle = 0° 




Fig. 1. — Diagram showing the definition of ^ used in the calculation of the forward-backward 
asymmetry for a single dec. band and a given 30 minute histogram. ^ is in the direction of 
hour angle. See text in Section 3.1 for details. 
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Fig. 2. — Representative plot of the number of events collected vs. azimuthal angle for 
the Milagro detector over the course of three days. The observed variation is the inherent 
geometric asymmetry in the detector which is at the level of ~ 10%. 
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Fig. 3. — Example histogram showing the number of cosmic-ray events vs. hour angle and 
declination containing seven years of data for a single sidereal half-hour. The solid circles 
are an example pair of pixels shown in Fig. 1 (with ^ = 42.5°) used in the calculation of the 
FB given by Eq.(l) & Eq.(3). 
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Fig. 4. — a) Sample 2-D histogram of FB vs. ^ and r.a. for a single dec. band (35° — 40°) 
calculated using Eq. [3] for seven years worth of data (the binning reflects that the data 
is collected over half hour intervals which are parameterized by ^o)- b) Sample histogram 
showing the result of the 2-D fit of Eq. IHto a) for a single slice in ^ = 40° — 45° (in black, left 
axis), and the reconstructed anisotropy of Eq. [2] using the six Fourier coefficients obtained 
by the 2-D fit to a) (in red, right axis). 
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Fig. 5. — Results of an all-dec-band anisotropy analysis using a MC generated UT signal. 
The signal consists of square well deficit from 150° to 210° in UT with an amplitude varying 
sinusoidally in time from 0.000 to 0.003. The signal seen in the ST and AST frames is 
induced by the UT time variation. The width of the curves reflects the statistical error. 
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Fig. 6. — Profiles in r.a. for individual 5° dec. bands from 10° to 50° used in the 2-D map 
seen in Figured The width of the lines reflects the statistical error. 
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Fig. 7. — Result of a harmonic fit to the fractional difference of the cosmic-ray rates from 
isotropic in equatorial coordinates as viewed by Milagro for the years 2000-2007. The color 
bin width is 1.0 x 10~^ reflecting the average statistical error. The two black lines show 
the position of the Galactic Equator and the solid circle shows the position of the Galactic 
North Pole. This map is constructed by combining 18 individual proflles of the anisotropy 
projection in r.a. of width 5° in dec. It is not a 2-dimensional map of the sky. The median 
energy of the events in this map is 6 TeV. 
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Fig. 8. — Amplitude and phase of the fundamental harmonic obtained from a fit to seven 
years of data for each 5° declination slice. The error bars are statistical. 




Fig. 9. — Anisotropy constructed by fitting a single projection containing data collected from 
all declinations collected over seven years in ST (red, left axis), UT (black, right axis) and 
AST (blue, right axis). The right axis (UT & AST) is expanded by a factor of 4 compared 
to the left axis (ST). The dashed curve is the predicted Compton-Getting (C-G) effect due 
to the Earth's motion around the sun which has an expected maximum value at 6 hr local 
solar time (196.3° UT given the longitude of Milagro). In these plots the width of the curve 
reflects the statistical error. 
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Fig. 10. — Total number of events over a six year period (left axis) and valley depth SBVD 
(right axis) vs. PMT multiplicity. The PMT multiplicity is defined as being the number of 
PMTs hit in the top layer of the pond for a given event. The lack of change in SBVD seen 
compared to the large decrease in the number of triggers shows explicitly the insensitivity 
of the analysis to variations in the trigger threshold. 
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Fig. 11. — Single-dec-band analysis showing seven years of data split into two sets: ^ < 25° 
(in yellow), and 25° < ^ < 50° (in green). Also plotted is the seven years of data split into 
three seasonal sets: Summer-Fall in black, Winter-Spring in blue, and Spring-Summer in 
red. The error bars are statistical + systematic. The independence of the analysis on ^ as 
well as the lack of seasonal variation can be seen here. 
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Fig. 12. — Position of minima in universal time vs. Modified Julian Date (MJD) for data 
collected over the 42 two-month periods starting in July 2000 and ending in July 2007. The 
error bars are the stat. errors. The solid line is the calculated position in universal time of 
a constant sidereal position at r.a. equal to 189°. The excellent agreement between the data 
and the calculation shows that the UT signal is dominated by the ST signal. 
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Fig. 13. — Histogram showing the projection of the AST curve from Figure [9] onto the 
anisotropy axis. The RMS of this distribution is 0.09 x 10"'^ which is used to determine the 
systematic error on the sidereal signaL 
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Fig. 14. — a) Sidereal time amplitudes of the three fit harmonics for the single band (all 
dec.) analysis, b) Sidereal time phases of the three fit harmonics for the single band analysis. 
Both plots contain seven yearly data sets from July 2000 to July 2007. The error bars are 
statistical. 
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Fig. 15. — Valley depth in the all-dec-band analysis (SBVD) vs. MJD for yearly sets from 
July 2000 to July 2007. The error bars are the hnear sum of the statistical & systematic 
errors. The solid line is the fit to a constant value and the dashed is the linear two-parameter 
fit. The x^/ndf for the fits are 86.2/6 and 4.4/5 respectively. The fit parameter in the fiat 
case is (2.39 ± 0.08) x lO"^; the two fit parameters to the function A{MJD) = po{MJD - 
53000) + pi are: po = (0.97 ± 0.11) x 10"^ and pi = (2.34 ± 0.08) x lO'^. 
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Fig. 16. — a) Universal time fit amplitudes for the single band (all dec.) analysis for seven 
yearly data sets from July 2000 to July 2007. b) Anti-sidereal time fit amplitudes for the 
single band analysis for yearly data sets. For the UT fundamental harmonic only we show 
the statistical error + an estimate of the systematic error. For AST the error bars are only 
statistical. Note the lack of any definite trend, as opposed to what is seen in ST (Fig. [14]) 
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Table 1. Fit parameters to the sidereal anisotropy for all 18 declination bands. 



Dec. Fund. Harm. 1st Harm. 2nd Harm. x^/d.o.f. Number 

(mean) Amplitude Phase Amplitude Phase Amplitude Phase of events 

(xlO-3) (deg) (xlO-3) (deg) (xlQ-^) (deg) (xlO^) 



77.5 


0.54±0.26 


2±27 


0.37±0.14 


13±11 


O.lOiO.ll 


-39±19 


262.57/234 


0.65 


72.5 


0.73±0.14 


22±11 


0.19±0.08 


-25±12 


0.06±0.06 


11±19 


266.50/282 


1.38 


67.5 


0.72±0.09 


23±7 


0.06±0.05 


-24±26 


0.01±0.04 


28±116 


308.67/330 


2.39 


62.5 


0.83±0.07 


19±5 


0.12±0.04 


-65±10 


0.15±0.04 


-2±4 


355.61/330 


3.63 


57.5 


0.99±0.06 


6±3 


0.12±0.03 


-42±8 


0.15±0.03 


3±4 


379.61/378 


4.98 


52.5 


1.10±0.05 


8±3 


0.22±0.03 


-60±4 


0.17±0.03 


6±3 


406.42/378 


6.31 


47.5 


1.31±0.04 


8±2 


0.33±0.03 


-63±2 


0.21±0.02 


2±2 


498.02/426 


7.51 


42.5 


1.71±0.04 


8±1 


0.44±0.02 


-68±2 


0.26±0.02 


1±2 


475.85/426 


8.46 


37.5 


1.95±0.04 


6±1 


0.45±0.02 


-73±1 


0.24±0.02 


3±2 


472.71/426 


9.07 


32.5 


2.04±0.04 


10±1 


0.47±0.02 


-76±1 


0.20±0.02 


1±2 


520.47/426 


9.28 


27.5 


2.17±0.04 


9±1 


0.53±0.02 


-78±1 


0.14±0.02 


0±3 


551.53/426 


9.07 


22.5 


2.39±0.04 


11±1 


0.52±0.02 


-81±1 


0.12±0.02 


-12±3 


564.14/426 


8.44 


17.5 


2.56±0.05 


12±1 


0.57±0.03 


-81±1 


0.10±0.02 


-32±5 


523.45/378 


7.45 


12.5 


2.62±0.06 


9±1 


0.61±0.03 


-85±2 


0.05±0.03 


-28±9 


397.37/330 


6.17 


7.5 


2.81±0.07 


5±1 


0.58±0.04 


-80±2 


0.08±0.03 


-39±7 


355.24/282 


4.74 


2.5 


3.05±0.10 


7±2 


0.61±0.05 


-80±3 


0.14±0.04 


54±5 


280.00/234 


3.31 


-2.5 


2.96±0.15 


8±3 


0.59±0.08 


-89±4 


0.23±0.06 


56±5 


271.14/186 


2.04 


-7.5 


3.80±0.38 


11±6 


0.42±0.19 


-85±13 


0.14±0.13 


52±18 


118.02/90 


1.04 



Note. — 



All quoted errors are statistical and are used in the calculation of x^. 
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Table 2. Fit parameters to the ST, UT and AST anisotropies using data collected over all 

declinations. 



Time 


Fund. Harm. 


1st Harm. 


2nd Harm. 




Frame 


Amplitude 


Phase 


Amplitude 


Phase 


Amplitude 


Phase 






(xlO-3) 


(deg) 


(xlO-3) 


(deg) 


(xlO-3) 


(deg) 




Sidereal 


1.994±0.012 


9.1±0.4 


0.400±0.007 


104.3±0.5 


0.118±0.006 


-0.9±1.0 


936.3/426 


Universal 


0.365±0.012 


217.4±1.9 


0.043±0.007 


168.7±4.8 


0.029±0.006 


93.2±4.3 


520.9/426 


Anti-Sidereal 


0.120±0.012 


47.1±5.8 


0.019±0.007 


155.5±11.1 


0.013±0.006 


4.9±9.4 


425.9/426 



Note. — All quoted errors are statistical and are used in the calculation of x^- For this single-dec-band 
analysis, the systematic errors exceed the statistical ones. 



